graph TD
A[Region 1]
A --> B[Clinic 1]
A --> C[Clinic 2]
A --> D[Clinic 3]
B -.-> M1[Patient 1]
B -.-> N2[Patient 2]
B -.-> O3[Patient 3]
C -.-> P4[Patient 4]
C -.-> Q5[Patient 5]
C -.-> R6[Patient 6]
D -.-> S7[Patient 7]
D -.-> T8[Patient 8]
D -.-> U9[Patient 9]
classDef region fill:#f9f,stroke:#333,stroke-width:2px;
classDef clinic fill:#ccf,stroke:#333,stroke-width:2px;
classDef patient fill:#cff,stroke:#333,stroke-width:1px;
class A,E,I region;
class B,C,D,F,G,H,J,K,L clinic;
class M1,N2,O3,P4,Q5,R6,S7,T8,U9,V10,W11,X12,Y13,Z14,A15,B16,C17,D18,E19,F20,G21,H22,I23,J24,K25,L26,M27 patient;
linkStyle default stroke-width:4px;
Mixed Models II
Eva Freyhult
NBIS, SciLifeLab
June 11, 2025
Nested and crossed designs
Grouping structures can be complex and consist of more than one level of grouping.
Nested: lower level units are only associated with one higher level unit
Crossed: lower level units are associated with multiple higher level units
Nested designs
In many cases, we have a hierarchical structure in our data, where observations are nested within groups. For example, patients within clinics and clinics within regions.
graph TD
E[Region 2]
E --> F[Clinic 4]
E --> G[Clinic 5]
E --> H[Clinic 6]
F -.-> V10[Patient 10]
F -.-> W11[Patient 11]
F -.-> X12[Patient 12]
G -.-> Y13[Patient 13]
G -.-> Z14[Patient 14]
G -.-> A15[Patient 15]
H -.-> B16[Patient 16]
H -.-> C17[Patient 17]
H -.-> D18[Patient 18]
I[Region 3]
I --> J[Clinic 7]
I --> K[Clinic 8]
I --> L[Clinic 9]
J -.-> E19[Patient 19]
J -.-> F20[Patient 20]
J -.-> G21[Patient 21]
K -.-> H22[Patient 22]
K -.-> I23[Patient 23]
K -.-> J24[Patient 24]
L -.-> K25[Patient 25]
L -.-> L26[Patient 26]
L -.-> M27[Patient 27]
classDef region fill:#f9f,stroke:#333,stroke-width:2px;
classDef clinic fill:#ccf,stroke:#333,stroke-width:2px;
classDef patient fill:#cff,stroke:#333,stroke-width:1px;
class A,E,I region;
class B,C,D,F,G,H,J,K,L clinic;
class M1,N2,O3,P4,Q5,R6,S7,T8,U9,V10,W11,X12,Y13,Z14,A15,B16,C17,D18,E19,F20,G21,H22,I23,J24,K25,L26,M27 patient;
linkStyle default stroke-width:4px;
Nested model
By using a nested model, we can account for the hierarchical structure of the data and identify the sources of variability in our data.
Nested models in R
In R and lme4, random effects for nested groupings can be included. Random intercepts for both regions and clinics nested within regions can be included as follows;
(1 | region/clinic) specifies that we have random intercepts for clinics nested within regions
Note, (1|region/clinic) is identical to writing (1|region) + (1|region:clinic).
A third way to implement the same model is to make the nesting implicit by coding the lower levels (here clinics) uniquely, so that a clinic name only exist in one region. With such a naming scheme, we can use the notation (1|region) + (1|clinic).
Nested models in R
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | region/clinic)
Data: data_int_nested
REML criterion at convergence: 1689.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.57686 -0.54281 -0.05204 0.44631 3.07737
Random effects:
Groups Name Variance Std.Dev.
clinic:region (Intercept) 473.98 21.771
region (Intercept) 2765.06 52.584
Residual 50.84 7.131
Number of obs: 240, groups: clinic:region, 12; region, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 105.7145 27.2460 3.880
x 0.3223 0.1617 1.993
Correlation of Fixed Effects:
(Intr)
x -0.124Crossed designs
In a crossed design, an observation is associated with multiple groups. For example, patients may be treated by different doctors at several clinics.
graph TD
A[Clinic 1]
E[Clinic 2]
I[Clinic 3]
B[Doctor 1]
C[Doctor 2]
D[Doctor 3]
M1[Patient 1]
M2[Patient 2]
M3[Patient 3]
M4[Patient 4]
M5[Patient 5]
A -.-> B
A -.-> C
A -.-> D
E -.-> B
E -.-> C
E -.-> D
I -.-> B
I -.-> C
I -.-> D
%% Patients connected to clinics to show they can come from any clinic, implying crossed with regions
B --> M1
B --> M2
B --> M3
B --> M4
B --> M5
C --> M1
C --> M2
C --> M3
C --> M4
C --> M5
D --> M1
D --> M2
D --> M3
D --> M4
D --> M5
%% Additional patient connections omitted for brevity
classDef region fill:#f9f,stroke:#333,stroke-width:2px;
classDef clinic fill:#ccf,stroke:#333,stroke-width:2px;
classDef patient fill:#cff,stroke:#333,stroke-width:1px;
class A,E,I region;
class B,C,D,F,G,H,J,K,L clinic;
class M1,M2,M3,M4,M5 patient;
linkStyle default stroke-width:4px;
Crossed models
A crossed design is coded like this in R:
y ~ x + (x | clinic) + (x | doctor)
Model selection in mixed models
Mixed models can become complex and model selection can be used to identify the best model.
Likelihood ratio tests (LRT) are often used to compare nested models, where one model is a simpler version of another. LRT can be used also for mixed models.
Mixed model have both random and fixed effects variables. We can use LRT to compare models with different fixed effects structures, but also different random effects structures.
What is LRT?
The likelihood is a measure of how well a model explains the data.
The likelihood ratio test compares the likelihoods of two models: the null model, which is a simpler model, and the alternative, more complex, model.
\[\lambda = - 2\ln{ \frac{L(restricted\,model)}{L(full\, model)}} \] The test statistic follows a \(\chi^2\)-distribution with degrees of freedom equal to the difference in the number of parameters between the two models.
A significant difference in likelihoods suggests that the more complex model explains the data better than the simpler model.
For mixed models, we can use the anova() function to perform LRT. We will start by testing the random effects structure, without change of fixed effects. Then we will test the fixed effects structure, keeping the random effects fixed.
ML and REML
In mixed models, we can fit mixed models using two different methods, Maximum Likelihood (ML) or Restricted Maximum Likelihood (REML).
ML has the advantage that it can be used to compare models with different fixed effects structures, but it gives biased estimates of the variance components. More about this in Payam’s maths section.
In summary, use
- ML when comparing models with different fixed effects structures (set
refit = TRUEinanova()) - REML when comparing diferent random effects in models models with the same fixed effects structure (set
refit = FALSEinanova()). Also, always use REML for the final model.
Random Effects Comparison
Sleep study data, random intercept.
sleepall <- sleepstudy |> filter(Days>=2) |> mutate(Days=Days-2)
m1 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepall)
summary(m1)Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
Data: sleepall
REML criterion at convergence: 1430
Scaled residuals:
Min 1Q Median 3Q Max
-3.6261 -0.4450 0.0474 0.5199 4.1378
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1746.9 41.80
Residual 913.1 30.22
Number of obs: 144, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 267.967 10.871 24.65
Days 11.435 1.099 10.40
Correlation of Fixed Effects:
(Intr)
Days -0.354Random Effects Comparison
Sleep study data, random intercept and random slope for Days.
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
Data: sleepall
REML criterion at convergence: 1404.1
Scaled residuals:
Min 1Q Median 3Q Max
-4.0157 -0.3541 0.0069 0.4681 5.0732
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 958.35 30.957
Days 45.78 6.766 0.18
Residual 651.60 25.526
Number of obs: 144, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 267.967 8.266 32.418
Days 11.435 1.845 6.197
Correlation of Fixed Effects:
(Intr)
Days -0.062Random Effects Comparison
Compare the two models using a likelihood ratio test.
Data: sleepall
Models:
m1: Reaction ~ Days + (1 | Subject)
m2: Reaction ~ Days + (1 + Days | Subject)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m1 4 1438.0 1449.9 -715.01 1430.0
m2 6 1416.1 1433.9 -702.05 1404.1 25.926 2 2.346e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fixed Effects Comparison
We can also compare models with different fixed effects structures. For example, we can compare a model with a fixed effect of Days to a model with a fixed effect of Days and Sex.
m2 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepall)
m3 <- lmer(Reaction ~ Days + sex + (1 + Days | Subject), data = sleepall)
anova(m2, m3, refit=TRUE)Data: sleepall
Models:
m2: Reaction ~ Days + (1 + Days | Subject)
m3: Reaction ~ Days + sex + (1 + Days | Subject)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m2 6 1425.2 1443.0 -706.58 1413.2
m3 7 1422.3 1443.1 -704.14 1408.3 4.8826 1 0.02713 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1P-value calculations
As you have noticed lme4 does not report p-values for the fixed effects.
In linear regression, p-values are calculated based on the t-distribution of the estimated coefficients. The calculation of denomicator degrees of freedom is straightforward since all observations are independent.
In mixed models, the situation is more complex due to the grouping structure. It is not straightforward to calculate the degrees of freedom for the fixed effects, since the observations are not independent.
There are however some approximation methods available to calculate the degrees of freedom and hence p-values for the fixed effects in mixed models.
We can use the lmerTest package, which extends lme4 to include p-values for fixed effects, based on e.g. Satterthwaite’s method.
lmerTest
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: conc ~ age + (1 | group)
Data: df
REML criterion at convergence: 304.1
Scaled residuals:
Min 1Q Median 3Q Max
-1.75481 -0.89181 0.09125 0.70848 1.85339
Random effects:
Groups Name Variance Std.Dev.
group (Intercept) 2224.44 47.164
Residual 85.47 9.245
Number of obs: 40, groups: group, 4
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 241.3334 26.0615 4.2275 9.260 0.00058 ***
age -1.8293 0.3014 36.1653 -6.068 5.54e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
age -0.422Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Reaction ~ Days + (1 + Days | Subject)
Data: sleepstudy
REML criterion at convergence: 1743.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.9536 -0.4634 0.0231 0.4634 5.1793
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 612.10 24.741
Days 35.07 5.922 0.07
Residual 654.94 25.592
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 251.405 6.825 17.000 36.838 < 2e-16 ***
Days 10.467 1.546 17.000 6.771 3.26e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
Days -0.138